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1.  Abstract 


This  study  presents  the  capability  of  the  use  of  the  peridynamic  laminate  theory  to 
capture  both  the  failure  progression  and  residual  strength  of  monolithic  and  composite 
laminates.  Predicting  damage  and  residual  strengths  of  composite  materials  involves 
capturing  complex,  distinct  and  progressive  failure  modes.  Peridynamics  is  a 
reformulation  of  classical  continuum  mechanics  that  utilizes  integral  equations  in  place  of 
partial  differential  equations  to  remove  the  difficulty  in  handling  discontinuities,  such  as 
cracks  or  interfaces,  within  a  body.  Damage  is  included  within  the  constitutive  model; 
initiation  and  propagation  can  occur  without  resorting  to  special  crack  growth  criteria 
necessary  in  other  commonly  utilized  approaches.  The  peridynamic  theory  realistically 
models  the  load  redistribution  arising  from  the  presence  of  complex  failure  modes 
through  the  use  of  multiple  interaction  types.  This  study  specifically  employs  an  inverse 
approach  to  obtain  the  critical  peridynamic  failure  parameters  necessary  to  capture  the 
residual  strength  of  a  structure.  The  validity  of  the  inverse  approach  is  demonstrated  by 
first  considering  its  application  in  determining  the  residual  strength  of  isotropic  materials 
with  pre-existing  cracks.  Its  validity  is  also  demonstrated  by  predicting  failure  loads  and 
final  failure  modes  in  a  laminate  with  various  hole  diameters  subjected  to  tensile  and 
compressive  loads. 

2.  Introduction 

By  utilizing  failure  analysis  simulations  that  are  able  to  capture  the  failure  modes  and 
structural  behavior,  the  residual  strength  of  structures  can  be  accurately  estimated  in  the 
design  phase.  Increased  confidence  in  these  simulations  enables  a  suitable  design  for  a 
structure  to  be  realized  more  efficiently  through  the  reduction  in  design  iterations  and 
number  of  experimental  test  required  for  structural  validation. 

The  ability  to  accurately  predict  the  failure  of  a  material  is  perhaps  one  of  the  most 
important  tasks  within  engineering,  as  understanding  the  failure  load  is  paramount  in 
designing  a  safe  structure.  Whereas  a  reasonable  estimate  can  be  safely  made  without 
understanding  the  complex  mechanisms  of  fracture,  lighter  structures  can  be  realized  with 
accurate  modeling  that  accounts  for  crack  initiation,  growth,  and  propagation  through  the 
media.  A  host  of  mechanisms  can  influence  the  initiation  of  failure,  as  well  as  the 
progression  through  a  structure.  Griffith  (1921)  investigated  failure  initiation  in  glass  by 
using  the  principle  of  conservation  of  energy,  which  helped  explain  the  large  discrepancy 
between  the  measured  stress  to  break  glass  and  the  theoretical  stress  required  to  break  the 
bonds  between  atoms.  This  became  the  basis  for  the  concept  of  Linear  Elastic  Fracture 
Mechanics  (LEFM).  Based  on  this  landmark  study,  numerous  investigations  were 
performed  to  explore  the  fracture  processes  and  establish  criteria  for  predicting  crack 
growth.  Williams  (1957)  showed  that,  within  the  realm  of  classical  continuum 
mechanics,  the  stress  field  near  the  crack  tip  approaches  infinity  in  an  elastic  and 
isotropic  material.  The  initiation  of  fracture  could  be  correlated  with  a  stress  intensity 
factor  to  describe  the  singular  stress  field.  The  presence  of  this  mathematical  singularity 
and  the  need  for  external  criteria  to  determine  failure  is  problematic  when  considering 
realistic  structures  with  complex  stress  fields. 
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Consideration  of  complex  geometries  and  loading  conditions  requires  the  interjection  of 
an  external  failure  criteria  into  a  suitable  framework  for  determining  how  and  when  a 
defect  in  the  form  of  a  crack  propagates.  The  finite  element  method  (FEM)  is  robust  in 
its  ability  to  determine  stress  fields  and  suitable  for  modeling  structures.  Different  failure 
criteria  were  incorporated  into  the  framework  of  FEM  with  varying  degrees  of  success.  It 
was  not  until  the  introduction  of  cohesive  zone  elements  by  Hillerborg  et  al.  (1976)  for 
Mode-I  fracture  mode  and  Xu  and  Needleman  (1994)  for  a  mixed-mode  fracture  that  a 
major  breakthrough  in  computational  fracture  mechanics  was  realized.  These  methods 
are  based  on  the  cohesive  zone  concept  independently  proposed  by  Dugdale  (1960)  and 
Barenblatt  (1962).  This  concept  introduces  a  restraining  stress  acting  across  the 
separating  surfaces  to  remove  the  stress  singularity  at  the  crack  tip.  Cohesive  zone 
elements  are  usually  surface  elements  that  are  placed  along  the  element  boundaries 
limiting  crack  growth  to  these  regions.  However,  the  crack  paths  are  unknown  a  priori, 
and  they  are  highly  sensitive  to  mesh  texture  and  alignment  (Klein  et  al.,  2001).  The 
concept  of  extended  Finite  Element  Method  (XFEM)  was  introduced  as  a  technique  to 
model  cracks  and  crack  growth  within  the  realm  of  finite  elements  without  remeshing 
(Belytschko  and  Black,  1999;  Moees  et  al.,  1999).  XFEM  is  based  on  the  partition  of 
unity  property  of  finite  elements  (Melenk  and  Babuska,  1996).  It  permits  cracks  to 
propagate  through  any  surface  within  an  element,  removing  the  limitations  of  the 
cohesive  zone  elements.  While  being  successfully  used  to  consider  numerous  fracture 
problems,  the  XFEM  still  requires  external  crack  growth  criteria  to  predict  crack  growth. 

When  considering  failure  in  complex  materials,  such  as  fiber  reinforced  composite 
laminates,  the  existing  failure  criteria  are  phenomenological  and  empirical  in  nature. 
There  exist  two  commonly  accepted  methods:  progressive  ply  failure  and  damage 
mechanics.  The  damage  mechanics  approach  employs  physically  based  equations  for 
damage  initiation  and  evolution  while  considering  the  microstructure  of  material  (Falzon 
and  Apruzzese,  2011;  Lapczyk  and  Hurtado,  2007;  Talreja,  1994).  However,  it  requires 
extensive  material  characterization  for  damage  parameters,  and  in  most  cases  such 
measurements  are  not  feasible.  Progressive  ply  failure  combines  failure  criteria  for 
damage  type  identification  (Hashin,  1980;  Sun,  2008)  and  degradation  of  the  material 
stiffness  (Chang  and  Chang,  1987;  Chang  and  Lessard,  1991;  Ochoa  and  Reddy,  1992). 
The  value  of  degradation  factor  is  not-physically  based;  it  is  assigned  a  small  enough 
value  in  order  to  ensure  the  convergence  of  finite  element  analysis  with  traditional  plate 
elements  for  modeling  composite  laminates. 

Silling  (2000)  and  Silling  et  al.  (2007)  introduced  the  peridynamic  theory  to  remove 
discontinuities  within  classical  continuum  theory  and  to  predict  failure  in  monolithic 
materials.  In  the  peridynamic  (PD)  theory,  internal  forces  are  expressed  through  nonlocal 
force  interactions  between  material  points  within  a  continuous  body.  Each  material  point 
interacts  with  other  material  points  within  a  finite  distance  referred  to  as  the  horizon,  and 
damage  is  part  of  the  constitutive  model.  The  resulting  equations  of  motion  do  not 
require  the  use  spatial  derivatives,  and  therefore,  are  defined  even  in  the  presence  of  a 
discontinuity. 
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The  PD  theory  is  a  reformulation  of  the  classical  continuum  equations  of  motion.  It 
replaces  the  partial  differential  equations  of  motion  with  integro-differential  equations. 
The  equations  of  motion  from  classical  continuum  theory,  can  be  expressed  as 

p(x)u(x,t)  =  V*o  +  b(x,t)  (1) 

in  which  p  is  the  mass  density,  and  o  is  the  first  Piola-Kirchoff  stress  tensor.  The 
displacement  and  body  force  density  vectors  are  denoted  by  u  and  b  respectively.  The 
position  vector  x  defines  each  point  with  respect  to  the  reference  configuration  at  time  t . 
In  the  PD  theory,  the  divergence  of  the  stress  tensor  on  the  right  hand  side  of  Eq.  (1)  is 
replaced  with  an  integral  that  accounts  for  all  the  nonlocal  force  interactions  that  exists 
between  a  material  point  and  all  of  the  material  points  within  its  horizon.  The  PD 
equations  of  motion  can  be  expressed  in  terms  of  the  force  density  vector,  t ,  that  exists 
between  material  points,  as 

/?(x)ii(x,t)  =  j*(t(u'-u,x'-x,t)-t'(u-u',x-x,,t))<i//+b(x,t)  .  (2) 

H 

where  the  material  point  x  interacts  with  material  point  x '  within  its  family  of  material 
points,  H .  Numerical  solution  to  the  resulting  equations  of  motion  can  be  achieved  by 
discretizing  in  the  form 


P(k)H(k)  1U(J)  U(k)’XU)  X(k)P) 

j= 1 

-KkXJ)  (u0->  -%)’XU)  -  w)]yo-)  +  bm  (3) 

With  the  inclusion  of  a  damage  parameter  in  the  force  density  vector  and  with  the 
removal  of  the  difficulties  related  to  discontinuities,  damage  initiates  and  progresses 
within  a  structure  where  it  is  energetically  favorable.  There  is  no  need  for  any  external 
failure  criteria  as  with  classical  continuum  based  methods.  PD  has  been  shown  to  be 
extremely  versatile  in  predicting  damage.  Silling  (2003)  considered  the  Kalthoff-Winkler 
experiment,  in  which  a  plate  having  two  parallel  notches  is  hit  by  an  impactor,  and 
showed  that  PD  is  able  to  predict  the  angle  of  crack  growth  observed  in  experiments. 
Impact  damage,  including  from  Charpy-V  notch  tests,  was  predicted  by  Silling  and 
Askari  (2004).  A  center  cracked  plate  was  used  by  Silling  and  Askari  (2005)  to  show 
numerical  convergence.  Gerstle  and  Sau  (2004)  demonstrated  the  ability  to  model 
damage  in  plain  and  reinforced  concrete  structures.  The  ability  to  incorporate  complex 
material  models  for  stretching  and  tearing  of  a  rubbery  material  was  demonstrated  by 
Silling  and  Bobaru  (2005).  Using  this  model,  they  were  able  to  predict  the  oscillatory 
path  of  a  blunt  tool  forced  through  a  membrane,  tearing  of  a  membrane,  and  the  crack 
growth  in  a  membrane  with  a  slit. 

In  recent  years,  PD  theory  has  also  been  applied  to  predict  damage  in  composites.  Within 
the  PD  framework,  the  simplest  approach  to  model  a  composite  layer  with  directional 
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properties  is  achieved  by  assigning  different  material  properties  in  the  fiber  and  other 
(remaining)  directions.  The  interactions  between  neighboring  layers  are  defined  by  using 
inter-layer  bonds.  Askari  et  al.  (2006)  and  Colavito  et  al.  (2007a,  2007b)  predicted 
damage  in  laminated  composites  subjected  to  low-velocity  impact  and  damage  in  woven 
composites  subjected  to  static  indentation.  In  addition,  Xu  et  al.  (2007)  considered 
notched  laminated  composites  under  biaxial  loads.  Also,  Oterkus  et  al.  (2010) 
demonstrated  that  PD  analysis  is  capable  of  capturing  bearing  and  shear-out  failure 
modes  in  bolted  composite  lap-joints.  Xu  et  al.  (2008)  analyzed  the  delamination  and 
matrix  damage  process  in  composite  laminates  due  to  low-velocity  impact.  Recently, 
Askari  et  al.  (2011)  considered  the  effect  of  both  high-  and  low-energy  hail  impacts 
against  a  toughened-epoxy,  intermediate-modulus,  carbon-fiber  composite.  Also,  Hu  et 
al.  (2011,  2012)  predicted  the  basic  failure  modes  of  fiber,  matrix,  and  delamination  in 
laminates  with  a  pre-existing  central  crack  under  tension.  The  analytical  derivation  of  the 
PD  material  parameters,  including  thermal  loading  conditions,  was  recently  given  by 
Oterkus  and  Madenci  (2012).  They  also  demonstrated  the  constraints  on  material 
constants  due  to  the  pair-wise  interaction  assumption.  An  alternative  approach  to  model 
composites  was  introduced  by  Kilic  et  al.  (2009)  by  distinguishing  fiber  and  matrix 
materials  based  on  the  volume  fraction.  Although  this  approach  may  bring  certain 
advantages  by  taking  into  account  the  inhomogeneous  structure  discretely,  it  is 
computationally  more  expensive  than  the  homogenized  approach.  Oterkus  et  al.  (2012) 
coupled  PD  with  FEM  to  predict  the  failure  loads  in  a  curved,  stiffened  composite  panel 
with  a  central  slit  subjected  to  uniaxial  loading  and  an  internal  pressure.  These  studies 
demonstrate  the  ability  of  the  PD  theory  to  accurately  model  both  the  progressive  damage 
and  final  failure  modes  of  composite  laminates. 

As  demonstrated  by  the  aforementioned  studies,  the  ability  of  PD  to  predict  failure 
initiation  and  progression  in  monolithic  and  composite  materials  has  been  well 
established.  However,  no  extensive  investigations  into  the  ability  of  PD  to  predict  the 
residual  strength  of  a  structure  exist.  Paramount  to  predicting  the  residual  strength  of  a 
structure  using  PD  is  the  assignment  of  critical  failure  parameters.  The  most  commonly 
used  failure  parameter  in  PD  is  the  critical  stretch.  The  stretch,  which  is  analogous  to 
strain  in  classical  continuum  mechanics,  is  monitored  between  PD  material  points.  When 
the  stretch  value  reaches  the  critical  stretch  then  the  interaction  between  those  material 
points  is  terminated.  Silling  and  Askari  (2005)  proposed  equating  the  energy  required  to 
create  a  unit  fracture  surface  in  PD  to  the  energy  release  rate  in  order  to  develop  an 
analytical  equation  to  determine  the  critical  stretch.  Alternatively,  Kilic  (2008)  showed 
that  for  monolithic  materials  that  the  failure  load  for  a  plate  with  a  central  hole  is  directly 
proportional  to  the  critical  stretch.  Based  on  these  two  concepts,  this  study  presents  an 
inverse  approach  to  determine  the  critical  stretch  value  required  for  the  accurate 
determination  of  the  residual  strength  of  a  structure.  It  applies  such  a  method  to  both 
monolithic  materials  and  fiber-reinforced  composite  laminates. 
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3 .  State-based  peridynamic  laminate  theory 

Although  the  use  of  laminated  composites  enables  the  design  and  fabrication  of  advanced 
airframes,  there  exist  significant  design  limitations  on  composite  structures.  These 
limitations  are  often  related  to  inter-laminar  stresses.  Thick  section  laminates,  especially 
in  the  transverse  direction,  experience  high  stresses  as  a  result  of  bending  moments, 
specific  local  stresses  and  impact  damage.  Many  applications  have  a  dynamic 
combination  of  high  load/high  cycle,  which  often  results  in  delamination  formation  and 
crack  growth.  Failure  of  composite  structures  involves  a  progressive  series  of  events 
with  discrete  failure  modes  such  as  matrix  cracking,  fiber-matrix  shear,  fiber  breakage, 
and  delamination.  Despite  the  development  of  many  important  concepts  to  predict 
material  behavior  and  failure,  the  prediction  of  failure  modes  and  residual  strength  of 
FRP  composites  is  still  a  challenge.  It  is  evident  that  the  inhomogeneous  nature  of 
composites  must  be  retained  in  the  analysis. 

The  peridynamic  theory  (PD)  pertaining  to  monolithic  materials  was  introduced  by 
Silling  (2000)  and  Silling  et  al  (2007).  Based  on  these  studies,  several  different 
approaches  to  model  composite  laminates  were  proposed.  In  the  PD  theory,  material 
points  interact  with  each  other  directly  through  the  prescribed  response  function,  which 
contains  all  of  the  constitutive  information  associated  with  the  material.  The  response 
function  includes  a  length  parameter  called  the  horizon,  8 .  The  locality  of  interactions 
depends  on  the  horizon,  and  interactions  become  more  local  with  a  decreasing  horizon. 
Each  fiber-reinforced  composite  lamina  of  a  laminate  shown  in  Fig.  1  is  idealized  as  a 
two-dimensional  structure  with  the  directional  dependency  of  the  interactions  between 
the  PD  material  points.  As  shown  in  Fig.  2,  the  material  point  ( q )  represents  material 
points  that  interact  with  material  point  ( k )  only  along  the  fiber  direction  with  an 
orientation  angle  of  6  in  reference  to  the  x-axis.  Similarly,  material  point  (r)  represents 
material  points  that  interact  with  material  point  ( k )  only  along  the  transverse  direction. 
However,  the  material  point  ( p )  represents  material  points  that  interact  with  material 
point  ( k )  in  any  direction,  including  the  fiber  and  transverse  directions.  The  orientation 
of  a  PD  interaction  between  the  material  point  ( k )  and  the  material  point  ( p)  is  defined 
by  the  angle  <f>  with  respect  to  the  x-axis.  The  domain  of  integration,  H  shown  in  Fig.  2 
is  a  disk  with  radius  8  and  thickness  h.  The  material  points  in  a  particular  lamina 
interact  with  the  other  material  points  of  immediate  neighboring  laminae  above  and 
below  it. 

As  shown  in  Fig.  1,  the  reference  coordinate  system  (x,y,z)  is  located  on  the  mid-plane 
of  the  laminate.  The  laminate  thickness,  h  is  given  by 

*  =  !>„  (4) 

n= 1 


where  N  is  the  total  number  of  lamina  in  the  stacking  sequence,  and  hn  is  the  thickness 
of  nth  lamina.  With  respect  to  the  mid-plane,  the  position  of  each  lamina,  z„  is  defined  as 
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(5) 


Figure  1.  Elevation  of  each  lamina  in  laminate  and  PD  material  points. 


Fiber  Direction,  F 
Transverse  Direction, T 
Arbitrary  Direction,  FT 


Figure  2.  PD  horizon  for  a  family  of  material  points  and  their  interactions  in  a  lamina. 


As  derived  by  Madenci  and  Oterkus  (2013),  the  equation  of  motion  for  material  point 
xn(k)  located  on  the  nth  layer  of  a  laminate  with  N  layers  can  be  expressed  as 


n  ••  ft  \  1  I  A_n  /  ,ft  ,,ft  it  _ft  ,  \  A_n  /,,ft  ,,ft  „ft  _ft  ,  \  I  t  /ft 

P(k)u(k)  ~  z-ftw)  (uo)  _u(i),x(j)  ~x(k)’t}~^(j)(k)  (U(i>  _u(j)’x(i)  ~xu),t)yu) 


7=1 


,  V  __  (ft)(m)T  7 m  .  ^  Y  n(nXmh/m  xh” 

Zj  F(^)  Zj  Z^Xi)  vOV 


(6) 


ftZ=ft+l  ,ft  — 


fti=ft+l,ft-l  7=1 
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where  the  material  point,  x"t)  on  the  n"‘  layer  is  associated  with  an  incremental  volume, 
V(nk)  and  a  mass  density  of  pn(k) ;  t  designates  time.  With  respect  to  a  Cartesian  coordinate 
system,  the  material  point,  x"k)  experiences  displacement,  u"t),  and  its  location  is 
described  by  the  position  vector,  y(k)  in  the  deformed  state.  The  displacement  and  body 
load  vectors  at  material  point,  x"t) ,  are  represented  by  u"t)  and  bn(k) ,  respectively.  The 
motion  of  a  material  point  conforms  to  the  Lagrangian  description  of  motion. 

Arising  from  in-plane  deformation,  tn(k)(j)  represents  the  force  density  that  material  point, 
x"  .}  exerts  up  on  material  point,  x"k  j .  The  force  density-stretch  relationships  can  be 
expressed  as 


U(*)’XG) 


y u)  y±k) 
y  u)  -y<*) 


(7) 


and 


t(y-X*)  K)  - u«)  •  *(*)  -  xO)  ‘  0  =  “  9  Bum  I  W_  0)| 

z  p»)  y (7>  i 

where  \k)(j)  and  B(j)(k)  are  auxiliary  parameters.  The  PD  strain  energy  density  at  point 
x(k)  for  in-plane  deformations  can  be  expressed  as 

w(k)=  a  e2(k)  - — - — r(|y0)  -  y (t)  I  -  |xu}  -  x(t)  |)2  vu) 

j= i  |x(y)  x(«| 

+  Kt  ^  lv  -v  |(ly(j)  “  y (t)  I "  lX(^>  “  I)' 

>=1  \x(j)  xm| 

+br  Zi — z — i(|yo)-y(«|-|x(;)-xw|)2yo)  (9> 

j=l  |X(y)  XG)| 

in  which  the  PD  parameter  a  is  associated  with  the  deformation  associated  with  the 
dilatation,  0(k),  of  the  matrix.  The  remaining  PD  parameters,  bF,  bT,  and  bFT,  are 

associated  with  deformation  of  material  points  in  the  fiber,  transverse,  and  arbitrary 
directions,  respectively.  The  fiber  and  transverse  interactions  exists  only  in  the  directions 
parallel  and  transverse  to  the  fiber,  and  therefore,  are  limited  to  j  =  \,J.  The  arbitrary 

interactions  exist  in  all  directions.  The  PD  dilation,  G(k) ,  of  the  matrix  can  be  expressed 
as 
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°o  ^ 

Qk)  3  I (1^0')  -yW|-|X0)  -X(t)D^-(tX;)^(y) 


i=1  XU)  X(*) 


The  auxiliary  parameters,  \k)(j)  and  Btj)(k) ,  can  be  determined  using  the  relationship 
between  force  density  vector  and  the  strain  energy  density,  W(k} ,  at  material  point  k , 
which  is  expressed  in  the  form 


*W0)(U0')  U(*)’ 


xo-)-x(«’ 


i  dw(k )  y  o-)-yw 


vo)5(|y0,-yw|)|y0)-yc 


Following  the  substitution  and  performing  the  differentiation,  the  auxiliary  parameters, 
\k)U)  and  Bl  j)(k) ,  can  be  expressed  in  terms  of  the  PD  material  parameters  as 


|  I  A(fcx./)^(fc)  ^^Bf^F  bpT  Bt^T  )  S(, 


(12a) 


Ix0)-xw| 


X«)~XU) 


r  A(jxk)fyj)  4b’ ( /-i/Bj,  +  bFT  +  fjjbj. )  y;)l 


(12b) 


v"  -v"  -  x"  -x" 

(«)  _K(y)  J(t)  A(i)  AW 

n  n 

XU)  ~X(k ) 


(13a) 


1  (x0)  -  x(Jt)  )//fiber  direction 


otherwise 


(13b) 


1  (x0)  -  x(k) )  _L  fiber  direction 


otherwise . 


(13c) 
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The  force  density  vectors,  and  with  m  =  (n  +  ]),(n-\)  develop  due  to  the 

transverse  normal  and  transverse  shear  deformations,  respectively,  between  the  material 
points  xn(k)  and  x"'( .  The  force  density-stretch  relationships  can  be  expressed  as 


0)0) 

V(k) 


k ) 


Vm  -v” 

*(*>  j(k) 

vm  -v” 

"(*)  J(*) 


(14a) 


and 


n( n)(m )  _ 


( n)(m )  y~(j)  -yw 

vm  -v" 

j(j)  "(*) 


(14b) 


where  C((^)(m) ,  and  D^u)  are  auxiliary  parameters.  Within  a  composite  laminate,  there 

commonly  exists  a  resin  rich  matrix  region  between  lamina  that  behaves  as  an  isotropic 
linear  elastic  material.  Therefore  the  strain  energy  expression  will  be  consistent  with  this 

form.  The  explicit  form  of  the  strain  energy  density  functions,  W"k)  and  W"k) ,  for 
transverse  normal  and  transverse  shear  deformations  can  be  written  as 


r  =h  V 

rr(k )  UN 


n=n+l,n—l  X(* ,  I 


yw-yw  -lx, 


(k) 


-y"  |Wm 

A(«|/  V(k) 


and 


<=ks  E  Zi 


m=n+l,n-l  j= 1  | X(j)  ^(k)  I 


-yo>  -ft 


w 


-0)1) 


vm 

v(j)  ’ 


(15a) 


(15b) 


in  which  the  PD  parameters  bN  and  bs  are  associated  with  the  transverse  normal  and 
transverse  shear  deformations,  respectively,  existing  between  bonded  lamina.  The 
parameter  8  is  the  horizon  size  in  the  thickness  direction  and  8  is  defined  as 

s=48r+¥. 


The  auxiliary  parameters,  ,  and  ,  can  be  determined  using  the  relationship 

between  force  density  vector  and  the  strain  energy  density,  W("k)  and  W"k) ,  at  material 
point  k  ,  which  are  expressed  in  the  form  (Madenci  and  Oterkus,  2013) 
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(16a) 


(n)(m) 

**(k) 


dWn 

uyv(k) 


y(*,-y; 


(*) 


y^-yw  yr*)-y"*) 


and 


q 


(k)(j) 


2  awfc  y”-)  -y"t) 

y(”  d(|yo)-ym|)|y^-y« 


(16b) 


Following  the  substitution  and  performing  the  differentiation,  the  auxiliary  parameters, 
C((^)(m) ,  and  can  be  expressed  in  terms  of  the  PD  material  parameters  as 


(n)(m)  _ 


y(») 


n  it  y (k)  y (k) 

<*>1/  m  _  » 

J(*)  "(*) 


(17a) 


and 


q 


( n)(jn ) 

(k)(j) 


=  4  b< 


1)1 

y o')  y (k) 

I'J 

yo-)-y« 

(17b) 


3.1  Material  parameters 

The  PD  material  parameters  can  be  obtained  in  terms  of  engineering  material  constants 
by  considering  simple  loading  conditions  and  equating  the  PD  strain  energy  density  to  the 
strain  energy  density  from  the  classical  continuum  mechanics.  The  PD  material 
parameters,  a,  d  ,  bF,  bT,  and  bFT ,  related  to  in-plane  deformations  can  be  obtained  by 

considering  four  different  loading  conditions:  simple  shear,  uniaxial  stretch  in  fiber 
direction,  uniaxial  stretch  in  transverse  direction,  and  biaxial  stretch.  The  PD  material 
parameters,  bN  and  bs ,  associated  with  transverse  deformations  can  be  obtained  by 

considering  two  different  loading  conditions:  transverse  normal  stretch  and  simple 
transverse  shear. 

As  derived  by  Madenci  and  Oterkus  (2013),  the  parameters  can  be  related  to  the  four 
independent  material  constants  of  elastic  modulus  in  the  fiber  direction,  En,  elastic 

modulus  in  the  transverse  direction,  in-plane  shear  modulus,  Gn,  and  in-plane 
Poisson’s  ratio,  v12  as 

«=^(er.-e«)  (isa) 
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(18b) 


d  = 


2 

7rhS2 


,  _  (2n  Qn  r^Q(,b) 

F~  (  .1 


28 

I 

\xn  -xn 

|A(7)  A(£) 

ly0) 

\J=l 

J 

(Q22  " 

2l2  2^66 ) 

(  N 

\ 

28 

a 

xn  -xn  \ 

A0')  A(£)| 

M'=1 

J 

b 

UFT 
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K+i+h, 


\3 


S2  +  2 


K+i+K 


s2  + 


K+l+K 
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+ 


I  2 


S2  +2 


K-i  +  K 
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K-\  +  K 
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where 


_  En  ^  _  vnE22  ^  _  E22  Q66  —  Gn 

^11  ’  ^12 


1  _  Vl2V2l 


1  -  V12V2\ 


1  “  ^12^21 


C  -_L  c  C  -_L  c  __L 

11  E  12  E  E  22  E  ’  66  G 

Ml  Ml  M2  M2  U12 


(18c) 

(18d) 

(18e) 

(18f) 


(18g) 

(19a) 

(19b) 


with  vnlEn=v2llE22  and  Em  and  Gm  represent  the  Young’s  modulus  and  shear 

modulus  of  the  matrix  material.  The  resulting  PD  model  of  a  laminate  consists  of 
laminae  connected  with  transverse  normal  and  shear  deformations,  and  it  accurately 
models  the  behavior  of  fiber  reinforced  laminate  composites. 
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The  PD  material  parameters  are  derived  under  the  assumption  that  the  material  point 
located  at  x  is  in  a  single  material  with  its  complete  neighborhood  entirely  embedded 
within  its  horizon,  5 .  However,  this  assumption  becomes  invalid  when  the  material 
point  is  close  to  free  surfaces  or  material  interfaces.  Therefore,  there  is  a  reduction  in 
material  stiffness  near  the  free  surfaces.  On  the  other  hand,  stiffness  near  the  interface  can 
exhibit  an  increase  or  reduction,  depending  upon  how  dissimilar  material  regions  interact 
across  their  interface.  Since  free  surfaces  and  material  interfaces  vary  from  one  problem 
to  another,  it  is  impractical  to  resolve  this  issue  analytically.  Therefore,  the  stiffness 
reduction  or  increase  due  to  surfaces  is  corrected  numerically.  A  surface  and  interface 
correction  factors,  as  explained  by  Madenci  and  Oterkus  (2013),  can  be  used  to  account 
for  the  loss  in  stiffness  near  the  surface  and  to  accommodate  the  transition  between 
dissimilar  laminae. 

A  numerical  treatment  is  necessary  to  solve  for  the  PD  equations  of  motion  which 
involves  differentiation  with  respect  to  time  and  spatial  integration.  The  numerical 
treatment  of  the  equations  of  motion,  Eq.  (6),  involves  the  discretization  of  the  domain  of 
interest  into  subdomains.  Within  these  subdomains,  the  velocity  and  displacement  fields 
are  assumed  to  be  constant.  Each  subdomain  can  then  be  represented  by  a  single 
collocation  point  located  at  the  center  of  each  subdomain.  The  details  of  the  numerical 
treatment  is  also  explained  by  Madenci  and  Oterkus  (2013). 

Solution  of  the  PD  equation  of  motion  requires  initial  velocity  and  displacement 
conditions  as  well  as  boundary  conditions.  The  introduction  of  displacement,  velocity,  or 
body  force  boundary  conditions  is  different  than  the  approach  used  in  classical  continuum 
mechanics  due  to  the  non-local  nature  of  PD.  The  details  of  the  introduction  of  the 
boundary  conditions  are  discussed  by  Madenci  and  Oterkus  (2013).  The  regions 
surrounding  the  boundary  conditions  require  special  treatment  to  prevent  unrealistic 
damage  accumulation.  No  fail  regions  are  introduced  in  the  region  surrounding  the 
boundary  conditions  to  prevent  this  phenomenon. 

Time  integration  of  Eq.  (6)  is  computationally  expensive  due  to  the  small  time-step  size 
required  for  stable  integration.  This  study  concerns  quasi-static  loading;  therefore,  the 
adaptive  dynamic  relaxation  (ADR)  method,  first  utilized  by  Kilic  and  Madenci  (2010),  is 
used  to  solve  the  system  of  ordinary  differential  equations  resulting  from  Eq.  (6). 


3.2.  Numerical  verification 

The  response  of  fiber-reinforced  composite  laminate  to  loading  is  dependent  on  the 
material  properties  of  each  lamina,  the  lamina  thicknesses,  and  the  specific  layup.  The 
validity  of  the  PDLT  for  composite  laminates  is  demonstrated  through  comparisons  with 
classical  laminate  theory  (CLT)  and  FEM  analysis  by  considering  laminates  with 
complex  layup  under  in-plane  tension  and  transverse  pressure. 
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3.2.1  Laminates  subjected  to  tensile  loading 

A  four  ply  laminate  subjected  to  uniaxial  tensile  loading  is  considered  to  demonstrate  the 
capability  of  PDLT  to  capture  the  deformation  of  a  laminate.  The  material  system  is 
fiberglass  epoxy  lamina.  The  material  properties  of  the  lamina  are  Eu  =38.6  GPa, 

E22  =  10.3  GPa  ,  Gu  =  4.5  GPa ,  and  v12  =  0.33 .  The  length  and  width  of  each  specimen 

is  L  =  1.6  m  and  W  =  1.0m,  respectively.  The  grid  size  in  the  model  is  Ax  =  0.01  mm 
resulting  in  160  points  in  the  length  direction  and  100  material  points  in  the  width 
direction.  The  nominal  ply  thickness  is  tk  =  0.01  m  resulting  in  a  total  laminate  thickness 
of  h  =  0.04  m .  The  horizon  is  specified  as  5  =  3.015Ax.  The  uniaxial  tension  is  applied 
as  a  body  load,  equivalent  to  a  stress  resultant  of  Nx  =1000N/m,  to  the  material  points 
in  a  volumetric  region  at  both  ends  of  the  plate  extending  the  width  of  the  laminate  and  a 
length  of  fc  =  0.1m.  During  the  solution,  failure  is  not  allowed  and  the  analysis  is  run 
until  the  system  reaches  equilibrium.  Four  separate  layups  are  considered  to  demonstrate 
the  ability  of  the  PD  formulation  to  capture  the  behavior  of  a  laminated  composite. 


b  b 


The  first  case  is  a  symmetric  cross-ply  laminate  with  a  layup  of  [0  /  90  ]s .  For 

symmetric  laminates,  CLT  predicts  that  there  is  no  coupling  between  bending  and 
extension.  Therefore,  the  panel  exhibits  only  in-plane  deformations  for  the  prescribed 
loading  condition.  Displacement  contours  on  the  mid-plane  of  the  laminate  are  shown  in 
Fig.  4.  The  variations  of  the  in-plane  displacement  components  at  the  mid-line  of  the 
laminate  are  compared  with  the  analytical  results  from  the  CLT  in  Figs.  5  and  6. 
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0.0 
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(a)  (b) 

Figure  4.  Displacement  contours  for  the  [0°  /  90°  ]s  laminate  subjected  to  an  uniaxial  tensile 

load:  a)  (, ux  )  and  b)  ( u ) 


Figure  5.  Horizontal  displacement  along  the  central  axis  for  the  [0°/90°]^  laminate. 
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Figure  6  Horizontal  displacement  (ii ,  )  perpendicular  to  the  central  axis  for  the  [0°  / 90  ]s 


laminate 


The  next  case  considers  a  unidirectional  [45°  ]4  angle  ply  laminate  with  the  fibers 

oriented  in  the  off-axis  direction  to  induce  a  stretching-shearing  coupling  behavior  when 
subjected  to  the  prescribed  loading  condition.  Figure  7  shows  the  stretching-shearing 
coupling  behavior  in  the  plane  of  the  laminate.  The  variations  of  the  in-plane 
displacement  components  at  the  mid-line  of  the  laminate  are  compared  with  the  analytical 
results  from  the  CLT  in  Figs.  8  and  9. 
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•  0.0 
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(a) 


(b) 


Figure  7  Displacement  contours  for  the  [45°  ]4  laminate  subjected  to  an  uniaxial  tensile  load:  a) 


(ux)  and  b)  (uy ) 
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Figure  8  Horizontal  displacement  ( ux  )  along  the  central  axis  for  the  [45°  ]4  laminate. 


Figure  9  Horizontal  displacement  (uy^j  perpendicular  to  the  central  axis  for  the  [45°  ]4 


laminate. 
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The  third  case  considers  an  anti-symmetric  [^45°  1-45°  /45°  /45°]  laminate.  When 
subjected  to  the  prescribed  uniaxial  tensile  loading,  it  experiences  stretching-twisting 
interaction  as  shown  in  Fig.  10.  The  variations  of  all  the  displacement  components  at  the 
mid-line  of  the  laminate  are  compared  with  the  analytical  results  from  the  CLT  in  Figs. 
11,  12,  and  13. 


(a)  (b)  (c) 

Figure  10  Displacement  contours  for  the  [—45  /  —45°  /  45°  /  45°]  laminate  subjected  to  an 
uniaxial  tensile  load:  a)  (ux ) ,  b)  (uy)  andc)  (//  ) 


ffi  -15 1- - 1 - 1 - 1 - 1 - *— 

-0.4  —0.2  0.0  0.2  0.4 

x— coordinate  (m) 


Figure  11  Horizontal  displacement  ( ux )  along  the  central  axis  for  the  [—45°  /  — 15  /  45°  /  45°] 

anti-symmetric  laminate. 
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Figure  12  Horizontal  displacement  [uy  )  perpendicular  to  the  central  axis  for  the 
[—45°  /  -4-5°  /  45°  /  45°  ]  anti-symmetric  laminate. 


Figure  13  Vertical  displacement  ( uz )  along  the  central  axis  for  the  [—45°  /  -45°  /  45°  /  45°  ] 

anti-symmetric  laminate. 
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The  final  case  considers  a  [0°  /  90°  /0°  /90°]  laminate,  which  exhibits  stretching-bending 
deformations  when  subjected  to  the  uniaxial  tension  as  shown  in  Fig.  14.  The  variation  of 
the  displacement  components  at  the  mid-line  of  the  laminate  are  compared  with  the 
analytical  results  from  the  CLT  in  Figs.  15,  16,  and  17. 

For  all  these  complex  laminate  layups,  the  PD  theory  captures  the  expected  deformation 
couplings,  and  there  is  a  good  agreement  between  the  PD  and  analytical  displacements. 


(a)  (b)  (c) 


Figure  14  Displacement  contours  for  the  [0°  /  90°  /  0°  /  90° ]  laminate  subjected  to  an  uniaxial 

tensile  load:  a)  (ux)  ,  b)  (//, )  andc)  (//, ) 


x— coordinate  (m) 


Figure  15  Horizontal  displacement  ( ux )  along  the  central  axis  for  the  [0°  /90°  /0°  /  90°  ] 

laminate. 
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Figure  16  Horizontal  displacement  [uy  )  perpendicular  to  the  central  axis  for  the 

[0°  /  90°  /  0°  /  90°  ]  laminate. 


Figure  17  Vertical  displacement  (uz )  along  the  central  axis  for  the  [0°  /  90°  /  0°  /  90° ] 

laminate. 
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3.2.2  Laminate  with  a  hole  subjected  to  tensile  loading 

A  laminate  with  a  layup  of  [0°  /  45°  /  45°  /0°] ,  with  a  central  hole,  is  subjected  to  uniaxial 
tensile  loading  to  demonstrate  the  capability  of  PD  laminate  theory  to  accurately  model 
complex  laminates.  The  material  system  is  unidirectional  AS4/3501-6  carbon  fiber 
epoxy  lamina.  The  material  properties  of  the  lamina  are  Eu  =  142  GPa ,  E22  =  10.3  GPa , 

G12  =  4.5  GPa ,  and  v12  =  0.27 .  The  length  and  width  of  each  specimen  is  L  =  1.0  m  and 

W  =  1.0  m,  respectively.  The  grid  size  in  the  model  is  Ax  =  0.01  mm  resulting  in  100 
points  in  the  length  direction  and  100  material  points  in  the  width  direction.  The  nominal 
ply  thickness  is  tk  =0.01  m  resulting  in  a  total  laminate  thickness  of  h  =  0.04  m.  The 
horizon  is  specified  as  c>=3.015Ax.  The  uniaxial  tension  is  applied  as  a  body  load, 
equivalent  to  a  stress  resultant  of  Nx  =4.26xl06  N/m,  to  the  material  points  in  a 
volumetric  region  at  both  ends  of  the  plate  extending  the  width  of  the  laminate  and  a 
length  of  b  =  0.03  m .  During  the  solution,  failure  is  not  allowed  and  the  analysis  is  run 
until  the  system  reaches  equilibrium. 

Figure  18  shows  the  contour  plots  of  the  displacement  components,  and  deformed  shape. 
The  variations  of  the  ux  displacement  in  the  x  direction  are  obtained  for  each  of  the 

layers  of  the  laminate  along  the  mid-line  as  well  as  along  the  edge  and  compared  with 
analytical  results  from  the  CLT  in  Figs.  19.  There  is  a  good  agreement  between  the  PD 
and  analytical  displacements  away  from  the  central  hole.  The  PDLT  clearly  captures  the 
coupling  between  in-plane  shear  and  extension  type  deformations. 


(a) 


(b) 


Figure  18  Displacement  contours  for  the  [0°  /  45°  /45°  /  0°]  laminate  subjected  to  an  uniaxial 

tensile  load:  a)  ( ux )  and  b)  {it  x ) 
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Figure  19  Horizontal  displacement  ( ux )  for  the  [0°  /  45°  /  45°  /0°]  laminate. 

3.2.3  Laminate  subjected  to  three  point  bending 

A  quasi-isotropic  laminate  with  layup  [0°  /90°  / -45°  /45°]5  is  subjected  to  three  point 

bending  to  demonstrate  the  capability  of  PDLT  in  the  presence  of  transverse  loading  as 
shown  in  Fig.20.  The  material  system  is  unidirectional  AS4/3501-6  carbon  fiber  epoxy 
lamina.  The  material  properties  of  the  lamina  are  En  =  142  GPa ,  E22  =  10.3  GPa  , 

Gn  =  4.5  GPa ,  and  v12  =  0.27  .  The  length  and  width  of  each  specimen  is  L  =  .1025  m 

and  W  =  .025  m,  respectively.  The  grid  size  in  the  model  is  Ax  =  0.0025  mm  resulting 
in  41  points  in  the  length  direction  and  10  material  points  in  the  width  direction.  The 
nominal  ply  thickness  is  tk  =  0.00013  m  resulting  in  a  total  laminate  thickness  of 
h  =  0.00104  m.  The  horizon  is  specified  as  5  =  3.0 1 5 Ax- .  A  displacement  boundary 
condition  of  u3=  0  was  enforced  on  the  material  points  residing  on  the  third  layer  along 
the  width  at  each  end  of  the  laminate.  The  central  load  is  applied  as  a  body  load, 
equivalent  to  a  pressure  of  Pz  =1000N/m2  ,  to  the  material  points  in  a  volumetric  region 
at  the  center  the  plate  extending  the  width  and  thickness  of  the  laminate  along  a  length  of 
b  =  0.02  m.  During  the  solution,  failure  is  not  allowed,  and  the  analysis  is  run  until  the 
system  reaches  equilibrium.  The  variations  of  the  uz  displacements  in  the  x  direction 

and  the  z  direction  are  obtained  at  the  mid-plane  of  the  laminate,  and  compared  with  the 
FEM  results  in  Figs.  21-22.  There  is  a  good  agreement  between  the  PD  predictions  and 
those  of  Refined  Zigzag  Theory  (RZT)  element  (Barut  et  al.  2013). 
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L  /2 


L  /2 


Figure  20  Loading  and  geometry  of  a  composite  laminate  under  three  point  bending. 


Figure  21  Vertical  displacement  (uz )  along  the  central  axis  for  the  [0°  /  90°  /  —45°  /  45  ]  v 

laminate. 
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Figure  22  Vertical  displacement  (uz )  along  the  vertical  axis  of  the  [0°  /  90°  /  —45°  /  45  ]s 

laminate. 


4.  Failure  prediction 

4.1  Damage  in  monolithic  material 

Material  damage  in  PD  is  introduced  through  elimination  of  interactions  among  the 
material  points.  It  is  assumed  that  when  the  stretch,  s(k)(j)  between  two  material  points  k 

and  j  exceeds  its  critical  value,  the  onset  of  damage  occurs.  Based  on  the  concept 

introduced  by  Silling  and  Askari  (2005),  the  critical  stretch  for  isotropic  materials  can  be 
obtained  as  derived  Madenci  and  Oterkus  (2013)  by  relating  the  energy  required  to  create 
a  fracture  surface  in  a  material  to  the  critical  energy  release  rate  of  that  material  as 
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It  is  worth  noting  that  the  critical  stretch  is  a  function  of  the  horizon.  The  value  of  the 
horizon  brings  in  the  effect  of  the  physical  material  characteristics,  nature  of  loading, 
length  scale,  and  the  computational  cut-off  radius.  This  simple  relationship  provides  the 
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value  of  critical  stretch  for  a  linear  elastic  brittle  material  with  a  known  critical  energy 
release  rate. 

In  order  to  include  damage  initiation  in  the  material  response,  the  force  density  vectors 
can  be  modified  through  a  history-dependent  scalar-valued  function  //  (Silling  and 
Bobaru,  2005).  When  the  stretch  between  these  material  points  exceeds  its  critical  stretch, 
failure  occurs;  thus,  the  history-dependent  scalar-valued  function,  //  is  zero  rendering  the 
associated  part  of  the  force  density  vector  to  be  zero.  Damage  is  reflected  in  the 
equations  of  motion  by  removing  the  force  density  vectors  between  the  material  points  in 
an  irreversible  manner.  As  a  result,  the  load  is  redistributed  among  the  material  points  in 
the  body,  leading  to  a  progressive  damage  growth  in  an  autonomous  fashion. 

Local  damage  at  a  point  is  defined  as  the  weighted  ratio  of  the  number  of  eliminated 
interactions  to  the  total  number  of  initial  interactions  of  a  material  point  with  its  family 
members.  The  local  damage  at  a  point  can  be  quantified  as  (Silling  and  Askari,  2005) 

|  /u(x'  -x,t)dV' 

— -  (21) 

H 

The  local  damage  ranges  from  zero  to  one.  When  the  local  damage  is  unity,  all  the 
interactions  initially  associated  with  the  point  are  eliminated,  while  a  local  damage  of 
zero  means  that  all  interactions  are  intact.  The  measure  of  local  damage  is  an  indicator  of 
possible  crack  formation  within  a  body.  For  example,  initially  a  material  point  interacts 
with  all  materials  in  its  horizon  as  shown  in  Fig.  23  (a);  thus,  the  local  damage  has  a 
value  of  0.  However,  the  creation  of  a  crack  terminates  half  of  the  interactions  within  its 
horizon  resulting  in  a  local  damage  value  of  Vi. 


(a)  (b) 

Figure  23  (a)  All  interactions  are  intact  (no  damage),  (b)  Half  of  the  terminated  interactions 

create  a  crack. 


At  each  time  step  in  the  numerical  process,  the  displacement  at  each  collocation  point  and 
the  stretch  between  collocation  points  is  computed.  The  stretch  is  monitored,  and  if  it 
exceeds  a  critical  value,  sc ,  then  the  interaction  is  terminated.  The  ratio  of  terminated 

interaction  to  total  interactions  at  each  collocation  point  is  calculated  during  the 
numerical  volume  integration  to  obtain  the  local  damage  defined  in  Eq.  (21). 
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A  fracture  surface  can  be  defined  as  a  region  where  no  interactions  between  material 
points  are  present.  Insertion  of  a  crack  in  the  PD  model  can  be  accomplished  by  breaking 
all  the  interactions  across  the  surface  of  the  crack.  The  insertion  of  a  crack  in  this  manner 
requires  determining  whether  an  interaction  crosses  the  plane  of  a  crack,  which  can  be 
tedious  depending  on  the  location  and  orientation.  If  the  interaction  crosses  the  crack 
surface,  then  the  interaction  is  terminated.  Alternatively,  a  crack  can  be  inserted  by 
terminating  the  interaction  if  one  of  the  material  points  is  above  the  plane  of  the  crack 
while  the  other  is  below  the  plane  of  the  crack. 

4.2  Damage  in  composite  laminates 

Five  different  interaction  types  are  present  in  the  PD  formulation  for  fiber  reinforced 
composite  laminates  to  capture  the  major  failure  modes.  This  is  accomplished  by 
assigning  unique  critical  stretch  values  to  the  different  interaction  types.  The  constitutive 
or  force-stretch  relations,  for  the  interactions  within  the  plane  of  a  lamina,  in  the  fiber, 
transverse,  and  arbitrary  directions,  are  shown  in  Fig.  24.  In  this  study,  the  transverse  and 
arbitrary  critical  parameters  are  combined  into  a  matrix  critical  stretch  in  tension  and 
compression  as  smt  and  smc ,  respectively.  The  critical  parameters  in  the  fiber  direction  in 

tension  and  compression  are  sft  and  sfc ,  respectively.  The  critical  stretch  parameters  for 

composite  laminates  can  be  obtained  using  experimental  methods  (Oterkus  et  al.,  2012), 
calibration  using  an  inverse  approach  (Colavito  et  ah,  2007a;  2007b)  or  through  equating 
the  energy  required  to  create  a  fracture  surface  to  the  energy  release  rate  (Hu  et  ah,  2011, 
2012). 


Interaction 


Figure  24  Force-stretch  relationships  for  peridynamic  interactions 

4.3  Inverse  approach  for  critical  stretch 

The  inverse  approach  is  an  alternative  to  the  use  of  analytical  critical  stretch  expression. 
This  method  combines  the  PD  simulation  with  an  experimental  test  to  calibrate  the 
critical  stretch  value.  For  the  first  step,  a  trial  critical  stretch  is  used  in  the  PD  simulation 
to  obtain  the  peak  failure  load  of  the  specimen.  The  PD  prediction  of  the  peak  failure 
load  is  then  compared  to  the  experimentally  observed  peak  failure  load  to  determine  if  the 
trial  critical  stretch  value  is  valid.  If  an  unacceptable  offset  exists  between  the  predicted 
and  measured  failure  loads,  then  the  trial  critical  stretch  is  adjusted  with  respect  to  the 
offset.  The  process  is  then  repeated  until  the  PD  prediction  of  failure  load  is  equal  to  the 
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experimental  failure  load.  At  the  end  of  the  iterative  process,  the  critical  stretch  value  is 
obtained. 

4.3.1  Monolithic  materials 

The  applicability  of  the  critical  stretch  as  a  failure  parameter  is  demonstrated  for  a  linear 
elastic  material  by  considering  the  experimental  study  conducted  by  Ayatollahi  and  Aliha 
(2009).  They  considered  diagonally  loaded  square  plate  specimens,  shown  in  Fig.  25a,  to 
investigate  the  effect  of  mode  mixity  ranging  from  pure  mode  I  to  pure  mode  II.  They 
provided  the  failure  loads,  crack  propagation  paths  for  each  of  the  specimens  and  fracture 
toughness  of  the  material,  KIC.  The  edge  length  of  the  diagonal  square  is  2 W  =  0.15  m 

and  its  thickness  is  h  =  0.005  m.  The  length  of  the  crack  is  2a  =  0.045  m  with  an 
orientation  angle  of  a .  The  material  has  an  elastic  modulus  of  E  =  2940  MPa , 

Poisson’s  ratio  of  v  =  0.38,  and  fracture  toughness  of  KIC  =  1.33  MPaVm  .  They  also 
reported  the  failure  loads  for  varying  crack  orientation  angles  of 
a  =  0°  (Mode I),  15°,  30°,  45°,  62.5° (Mode II)  .  Center  of  the  crack  coincides  with  the  origin 
of  the  Cartesian  coordinate  system.  The  applied  load  is  introduced  through  a  velocity 
constraint  of  V7  =  1  x  1  IT9  m/s  in  the  volume  defined  by  dotted  lines  while  enclosing  the 
circular  regions  at  opposite  corners.  The  initial  crack  is  inserted  in  the  PD  model  by 
removing  the  interactions  across  the  crack  surface. 


Figure  25  a)  DLSP  Specimen  b)  PD  model  of  DSPL. 

The  crack  propagation  paths  obtained  from  the  PD  simulations  and  those  of  the 
experimental  observations  are  plotted  in  Fig.  26.  A  good  agreement  is  observed  between 
the  predictions  and  observed  crack  paths.  Crack  growth  initiation  angles  are  also 
compared  between  the  predictions  and  measurements.  Again,  a  very  good  comparison  is 
obtained  as  shown  in  Fig.  27.  Finally,  the  failure  loads  are  compared,  and  it  is  observed 
that  the  failure  loads  obtained  from  the  PD  simulations  are  within  15%  of  the 
experimental  values  for  all  crack  inclination  angles  as  depicted  in  Fig.  28.  While  the  PD 
simulations  are  close  to  the  experimental  results  for  pure  Mode  I  case  and  pure  Mode  II 
cases,  the  mixed  mode  PD  failure  predictions  are  higher  than  the  experimental  values.  A 
possible  reason  for  this  could  be  due  to  specimen  preparation,  which  does  not  ensure  a 
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sharp  crack  tip.  The  specimen  preparation  procedure  requires  the  use  of  a  fret  saw  to 
create  a  crack  without  additional  cyclic  loading.  The  inclined  angle  of  the  crack,  coupled 
with  the  shape  of  the  crack  tip,  could  in  effect  change  the  cracks  tip  orientation  causing 
the  offset  observed  in  the  comparison.  Despite  this  offset,  the  agreement  is  satisfactory; 
thus,  validating  the  critical  stretch  value  of  0.089. 


Figure  26  Experimental  and  PD  crack  propagation  paths 
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Figure  27  Comparison  of  crack  growth  initiation  angle  between  peridynamic  and  experimental 
results  as  a  function  of  crack  inclination  angle. 


30 


3500 


3000 


2500  i 


"§  2000 k 

_o 

3  15001- 

•<s> 

£ 

1000  - 


500  - 


A  ▲  Experimental  values 
•  •  PD  predictions 


15  30  45 

Crack  inclination  angle  ( degrees ) 


62.5 


Figure  28  Comparison  of  the  failure  load  between  peridynamic  and  experimental  results  as  a 

function  of  crack  inclination  angle. 

4.3.2  Laminated  composites 

The  inverse  approach  for  monolithic  materials  is  modified  for  the  PDLT  to  accommodate 
the  multiple  interaction  types  and  critical  stretch  values.  To  demonstrate  the  approach, 
the  inverse  method  is  used  to  obtain  the  fiber  and  matrix  critical  stretch  values  for  a 
notched  quasi-isotropic  AS4/3501-6  laminates  subjected  to  tensile  and  compressive 
loading.  In  the  experimental  study  conducted  by  Wang  et  al.  (2004),  the  tensile  and 
compressive  peak  loads  are  obtained  for  a  [45  /  0/-45/  90] ls  panel  with  four  different 

central  hole  diameters  (2.00,  3.81,  6.35,  and  9.55  mm).  The  following  procedure  is  used 
to  obtain  the  matrix  and  fiber  critical  stretch  values  for  a  horizon  value  of  0.0022914  m. 

To  utilize  the  inverse  approach,  several  assumptions  regarding  the  progression  of  failure 
are  necessary.  The  first  assumption  is  that  the  failure  in  tension  depends  primarily  on 
matrix  cracking.  Therefore,  the  peak  failure  load  of  the  laminate  primarily  depends  on 
the  tensile  matrix  critical  stretch.  This  assumption  is  made  after  examination  the  failure 
patterns  for  two  different  laminates,  [0/90]^  and  [45/ -45]^.  The  failure  modes 

present  in  the  [0  /  90]  ^  laminate  are  longitudinal  fracture  of  the  matrix  in  the  90°  layer 
and  fracture  of  the  fibers  in  the  0°  layer.  The  failure  modes  present  in  the  [45  /  —45]^ 

laminate  are  longitudinal  fracture  of  the  matrix  for  all  layers  and  delamination  between 
layers.  For  both  laminates,  the  longitudinal  fracture  of  the  matrix  precedes  the  other 
failure  modes  and  dictates  the  final  failure  pattern,  and  allows  for  such  an  assumption. 
The  second  assumption  concerning  the  laminate  compressive  strength  is  dictated  by 
compressive  fiber  critical  stretch.  This  assumption  relies  on  the  experimental 
observations  of  the  progressive  failure  of  the  identical  laminate  specimen  with  a  6.35  mm 
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diameter  hole.  Suemasu  et  al.  (2005)  observed  that  damage  first  appeared  due  to  fiber 
micro-buckling  in  the  0°  layers.  They  also  observed  additional  fiber  micro-buckling  in 
the  0°  layers  and  inter-laminar  delamination  prior  to  a  sudden  failure  of  the  laminate. 
Finally,  the  relationship  between  compressive  and  tensile  critical  stretch  values  is 
assumed  to  be  equivalent  to  the  ratio  of  the  compressive  and  tensile  strengths  of  the 
lamina.  The  longitudinal  and  compressive  strength  ratios  are  assumed  to  be  equivalent  to 
the  fiber  and  matrix  critical  stretch  ratios,  respectively.  The  first  two  assumptions  allow 
the  tensile  matrix  and  compressive  fiber  critical  stretch  values,  while  the  third  assumption 
allows  the  compressive  matrix  and  tensile  fiber  critical  stretch  values  to  be  obtained. 

The  laminate  with  the  6.35  mm  central  hole  is  used  to  calibrate  the  critical  stretch  values 
due  to  the  additional  experimental  observations.  In  the  first  step,  the  tensile  matrix 
critical  stretch  is  calibrated  using  iterative  PD  simulations  of  the  tensile  loading  of  the 
notched  laminate.  The  tensile  matrix  critical  stretch  is  adjusted  until  the  peak  failure  load 
for  the  PD  simulation  matches  the  measured  value.  For  this  step,  the  tensile  fiber  critical 
stretch  is  assigned  an  artificially  high  value  to  prevent  failure  due  to  fiber  damage. 
Following  this  calibration,  the  tensile  matrix  critical  stretch  value  is  used  to  obtain  the 
compressive  matrix  critical  stretch  value  using  the  ratio  of  tensile  to  compressive 
transverse  strength  of  the  lamina,  Yt  and  Yc,  respectfully.  The  compressive  matrix 

critical  stretch  is  obtained  by  multiplying  the  strength  ratio  (Yc/Yt)  by  the  matrix  tensile 

critical  stretch.  Then,  the  compressive  fiber  critical  stretch  value  is  obtained  in  a  similar 
manner  using  iterative  PD  simulations  of  the  notched  laminated  subjected  to  compressive 
loading  and  adjusting  the  compressive  fiber  critical  stretch  to  obtain  equivalent  failure 
loads.  Then,  the  calibrated  compressive  fiber  critical  stretch  value  can  be  used  to  obtain 
the  tensile  fiber  critical  stretch  value  using  the  ratio  of  tensile  to  compressive  longitudinal 
strength  of  the  lamina,  Xt  and  Xc ,  respectfully.  The  entire  process  is  repeated  to  fine- 

tune  the  calibration  of  the  critical  stretch  values  until  the  failure  loads  approach  the 
measured  values. 

4.3.2. 1  Laminated  plate  with  a  center  hole 

The  quasi-isotropic  AS4/3501-6  carbon/epoxy  composite  laminates  with  a  hole  were 
previously  considered  by  Wang  et  al.  (2004).  The  layup  for  all  of  the  composite  panels  is 
[45/0/^45  /  90]2S .  Each  laminate  specimen  is  305  mm  x  38.1  mm  x  2.08  mm  with  a 

central  through  thickness  hole.  The  specimen  dimensions  are  shown  in  Fig.  29.  Four 
different  diameters  ( 2.00,  3.81,  6.35,  and  9.55  mm)  of  central  holes  are  considered.  The 

material  properties  of  the  lamina  are  Eu  =  142  GPa ,  E22  =  10.3  GPa ,  Gu  =  4.5  GPa ,  and 
v12  =  0.27 .  The  strength  properties  of  the  lamina  are  erlt  =  2280  MPa ,  <rlc  =  -1440  MPa 
,  <J2t  =  57  MPa ,  a2c  =  -228  MPa ,  and  r12  =  100  MPa.  The  laminated  specimens  were 

tested  to  failure  under  both  tensile  and  compressive  loading  in  accordance  with  ASTM 
3039  and  SACMA  SRM  3R  standards,  respectively. 
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Experimental  Specimen: 


305mm 

(2.0mm,  3.81mm,  _ 

40 

6.3  5mm,  &  9.55mm) 

38.1mm 

< - 

Peridynamic  Model: 


118.5mm 

(2.0mm,3.81mm,  _ 
6.3  5mm,  &  9.55mm) 

-0 

38.1mm 

Figure  29.  Dimension  and  loading  conditions  for  notched  quasi-isotropic  AS4/3501-6  laminates. 

The  dimensions  of  each  PD  model  are  1 18.5  mm  x  38.1  mm  with  a  total  thickness  of  2.08 
mm  as  shown  in  Fig.  29.  Each  lamina  is  modeled  with  a  grid  of  156  material  points  in 
the  length  direction  and  50  material  points  in  the  thickness  direction;  the  grid  spacing  in 
the  plane  of  the  lamina  is  0.76  mm.  The  horizon  radius  is  specified  as  8  =  3.015Ax.  A 
central  hole  is  added  to  each  lamina  by  removing  the  material  points  falling  within  its 
specified  diameter.  The  spacing  in  the  thickness  direction  between  each  lamina  is  0.13 
mm.  The  loading  is  applied  by  specifying  an  end  displacement  at  both  the  top  and 
bottom  edge  of  the  model  in  incremental  steps  of  0.5  mm.  Convergence  for  each  step  is 
ensured  to  enforce  the  quasi-static  loading  conditions. 

The  critical  stretch  parameters  are  obtained  through  the  inverse  approach  by  calibrating 
the  experimentally  observed  peak  failure  loads  to  those  obtained  from  the  PD  simulation. 
This  calibration  was  completed  on  the  specimen  with  a  6.35  mm  diameter  central  hole. 
The  tensile  matrix  critical  stretch  and  compressive  fiber  critical  stretch  values  are 
obtained  through  calibration  with  the  tensile  and  compressive  failure  loads,  respectively. 
The  ratios  of  compressive  strengths  for  a  lamina  are  then  used  to  determine  the  other 
critical  stretch  values.  The  critical  stretch  values  obtained  using  the  inverse  approach  are 
scmt  =  0.0176,  scmc  -  -0.0528,  ^=0.01882,  and  scfc  =  -0.01 1 89 .  A  PD  simulation  of  a 

laminate  with  a  6.35  mm  diameter  central  hole  subjected  to  tensile  loading  is  used  to  test 
the  ability  of  these  values  to  capture  the  failure  load.  The  simulation  predicted  the 
initiation  of  failure  at  the  central  hole  followed  closely  by  an  accumulation  of  damage 
throughout  the  entire  laminate. 
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The  peak  loads  obtained  from  the  PD  simulations  are  within  12%  of  the  experimentally 
observed  failure  loads.  Figure  30  shows  the  peak  failure  loads  obtained  from  the  PD 
simulations  along  with  the  experimentally  results  reported  by  Wang  et  al.  (2004)  for  the 
four  different  hole  diameters  that  were  considered.  A  good  agreement  is  observed 
between  the  experimental  and  the  PD  predicted  peak  failure  loads  for  both  the  tensile  and 
compressive  tests. 


(a)  (b) 

Figure  30.  Experimental  (Wang  et  al.,  2004)  and  PD  peak  loads  for  a)  tensile  and  b)  compressive 

loading. 


Tensile  Loading  -  The  final  damage  pattern  transitions  from  one  that  is  aligned  with  the 
45°  angle  with  respect  to  the  principle  coordinate  system  for  the  small  diameter  hole  size, 
to  the  one  aligned  with  the  90°  angle  for  the  larger  hole  sizes.  Figure  31  shows  the  final 
damage  patterns  for  each  ply  in  the  laminate  composite  subjected  to  a  tensile  load  for 
four  different  central  hole  diameters.  The  specimen  fracture  pattern  can  be  obtained  by 
plotting  the  minimum  damage  at  each  material  point  in  the  thickness  direction.  The 
fracture  for  each  of  the  specimens  is  plotted  in  Fig.  32.  The  resulting  fracture  pattern 
transition  from  a  slant  to  a  flat  failure  mode  as  the  hole  size  increases.  A  similar 
observation  shown  in  Fig.  33  was  captured  by  Poon  (1991)  for  carbon  fiber  epoxy 
specimens  with  the  same  base  layup  pattern  and  three  times  the  thickness. 

The  initial  tensile  failure  mode  in  all  of  the  specimens  is  matrix  cracking  originating  at 
both  sides  of  the  central  hole  in  the  90°  plies  closest  to  the  surface  of  the  laminate.  As 
the  loading  increases,  the  damage  accumulates  in  all  of  the  layers  and  propagates  towards 
the  edges  of  the  plate.  Delamination  due  to  shear  is  observed  between  the  layers  near  the 
surface.  The  amount  of  matrix  and  delamination  that  occurs  depends  on  the  size  of  the 
central  hole.  In  the  specimens  with  a  large  central  hole,  the  amount  of  damage  that 
accumulates  is  focused  on  either  side  of  the  hole.  In  the  specimens  with  a  small  central 
hole,  the  area  where  damage  occurs  is  extensive.  The  fiber  failure  occurs  first  in  the  0° 
plies,  followed  quickly  by  the  +  45°  layers.  The  peak  load  is  observed  with  the  initiation 
of  the  fiber  failure. 
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Figure  31..  Failure  patterns  in  each  ply  of  a  [45  /  0  /  —45  /  90]2^  laminate  with  a  central  hole  of 
diameter  a)  2.00  mm,  b)  3.81  mm,  c)  6.35  mm,  and  d)  9.55  mm  subjected  to  tensile  loading. 


Figure  32.  Failure  patterns  for  a  [45  /  0  /  —45  /  90]2S  laminate  with  a  central  hole  of  diameter  a) 
2.00  mm,  b)  3.81  mm,  c)  6.35  mm,  and  d)  9.55  mm  subjected  to  tensile  loading. 
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Figure  33  Effect  on  failure  patterns  for  a  [45  /  0  /  —45  /  90]65  laminate  with  a  central  hole  of 

due  to  increase  in  hole  size  (Poon  1991). 


Compressive  Failure  -  Figure  34  shows  the  final  damage  patterns  for  each  ply  in  the 
laminate  composite  subjected  to  a  compressive  load  for  the  four  different  diameter  central 
hole  sizes  considered.  Due  to  quick  propagation  in  the  cases  with  the  small  diameter 
holes  (2.00  mm  and  3.81  mm),  the  final  damage  pattern  is  obtained  at  a  time  step  long 
after  the  peak  failure  load.  In  all  cases,  the  damage  extends  horizontally  from  the  central 
hole  toward  the  edge  of  the  specimen.  The  final  fracture  pattern  at  each  material  point  in 
the  thickness  direction  is  plotted  in  Fig.  35  for  each  of  the  different  hole  diameters.  The 
fracture  pattern  is  consistent  for  all  the  different  hole  diameters. 
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Figure  34  Failure  patterns  in  each  ply  of  a  [45  /  0/  —45  /  90]25  laminate  with  a  central  hole  of 
diameter  a)  2.00  mm,  b)  3.81  mm,  c)  6.35  mm,  and  d)  9.55  mm  subjected  to  compressive  loading. 


Figure  35  Failure  patterns  for  a  [45  /  0  /  —45  /  90]25  laminate  with  a  central  hole  of  diameter  a) 
2.00  mm,  b)  3.81  mm,  c)  6.35  mm,  and  d)  9.55  mm  subjected  to  compressive  loading. 
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The  PD  simulation  for  progressive  failure  can  be  compared  with  the  experimental 
observations  (Fig.  36)  obtained  by  Suemasu  et  al.  (2006).  The  initial  damage  observed  in 
the  PD  simulation  is  fiber  failure  due  to  compression  in  the  0°  layers  followed  by  shear 
failure  leading  resulting  in  the  loss  of  stiffness  of  the  structure.  This  corresponds  with  the 
observations  by  Suemasu  et  al.  (2006)  of  an  the  initial  failure  mode  consisting  of  fiber 
micro-buckling  in  the  0°  layers  followed  by  delamination  immediately  before  the  failure 
of  the  specimen.  The  PD  simulation  also  captures  the  damage  extending  horizontally 
from  the  hole  that  is  observed  in  a  C-scan  of  the  experimental  specimen  immediately 
before  the  final  failure  of  the  specimen. 


Figure  36  a)  C-Scan  of  the  damage  immediately  before  final  failure  (Suemasu  et  al.,  2006),  b) 
Failure  pattern  predicted  using  PD  simulation. 


5.  Final  remarks 

This  study  demonstrates  the  ability  of  the  PDLT  to  accurately  capture  the  material 
response  and  progressive  failure  in  composite  laminates.  PDLT  correctly  captures  the 
structural  response  of  laminates  under  various  loading  conditions  and  laminate  layups. 
The  capability  of  the  PDLT  to  accurately  model  complex  laminates  is  verified  against  the 
classical  laminate  theory  and  FEM. 

The  PDLT  is  able  to  capture  matrix  cracking,  fiber  breakage,  and  delamination  due  to 
transverse  normal  and  transverse  shear  deformation.  Also,  the  PD  simulations  using  the 
critical  stretch  values  obtained  using  the  inverse  approach  are  able  to  capture  the  residual 
strength  of  the  notched  specimens  subjected  to  both  tensile  and  compressive  loading. 
The  use  of  PD  in  this  methodology  allows  the  progressive  damage  of  the  composite 
structure  to  be  captured  in  an  extremely  accurate  manner. 
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